lisa.perm<-function(x,listw,perm=999){
    
    z<-x-mean(x)
    n<-length(x)
    s2<-sum(z^2)/(n-1)
    lz<-lag.listw(listw,z)
    Ii<-(z/s2)*lz
    
    lisas<-matrix(0,n,perm)
    rids<-matrix(0,max(card(listw$neighbours)),perm)
    n_1<-n-1
    
    for(j in 1:perm){
        rids[,j]<-sample.int(n_1,size=max(card(listw$neighbours)),replace=FALSE)
    }
    
    for(i in 1:n){
        #	x=sample(z[-i],replace=FALSE) 
        x<-z[-i]
        rz<-x[as.numeric(rids[1:length(listw$weights[[i]]),])]
        rz<-matrix(rz,nrow=length(listw$weights[[i]]),ncol=perm,byrow=FALSE)
        lisas[i,]=z[i]*(listw$weights[[i]]%*%rz)
    }
    
    rlisas<-(1/s2)*lisas
    
    above<-rlisas>=Ii
    larger<-rowSums(above)
    low.extreme<-(perm-larger)<larger
    larger[low.extreme]<-perm-larger[low.extreme]
    
    p.sim<-(larger+1)/(perm+1)
    
    y<-cbind(Ii,p.sim,z)
    colnames(y)<-c('Ii','p-val','Zi')
    class(y)<-'lisaPerm'
    return(y)
}


plot.lisaPerm<-function(lmoran,shape,signif=.05,legend.title=''){
    if(class(lmoran)!='lisaPerm')stop('object class is not lisa.perm')
    
    quadrant <- vector(mode="numeric",length=nrow(lmoran))
    
    quadrant[lmoran[,3] >0 & lmoran[,1]>0] <- 1 #H-H   
    quadrant[lmoran[,3] <0 & lmoran[,1]>0] <- 2 #L-L
    quadrant[lmoran[,3] <0 & lmoran[,1]<0] <- 3 #L-H     
    quadrant[lmoran[,3] >0 & lmoran[,1]<0] <- 4 #H-L
    quadrant[lmoran[,2]>signif] <- 0    
    
    brks <- c(0,1,2,3,4)
    palette <- c("white","red","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4))
    colors<-palette[findInterval(quadrant,brks,all.inside=FALSE)]
    
    plot(shape@bbox,t='n',axes=FALSE,ylab='',xlab='')
    for(i in 1:nrow(shape)){
        polygon(shape@polygons[[i]]@Polygons[[1]]@coords,col=colors[i],lty=3)
    }
    box()
    legend("bottomright",legend=c("insignificant","high-high","low-low","low-high","high-low"),
           fill=palette,bty="n",cex=0.9,y.intersp=1,x.intersp=1,title=legend.title)
    title("LISA Cluster Map")
    
}